Lesson 4 - Interpret Output

Calibration Diagnostic:

For each iteration of the calibration all the metrics listed in Table below are calculated and a number of plots are generated to ease the assessment of the progress made at a given iteration. These plots help to check whether the calibration is heading toward the right direction and diagnose for possible bugs and disruption in the calibration process. The table below summarized all the metrics that are currently calculated as a part of calibration process.

Metric Equation Optimal Value Purpose
Objective Function User specified 0 Used in the minimization process of DDS
Nash-Sutcliffe Efficiency (NSE) $ NSE = 1 - \frac{\sum (Q_o - Q_s)^2 } {\sum (Q_o - \bar{Q_o})^2} $ 1 Single metric combining timing and magnitude errors
Normalized Nash-Sutcliffe Efficiency (NNSE) $ NNSE = \frac{1}{2-NSE} $ 1 Normalized Single metric combining timing and magnitude errors,
Normalized Nash-Sutcliffe Efficiency Squared (NNSEsq) $ NNSEsq = \frac{1}{2-NSE^{Sq}} \\ where: \\ NSE \space streamflow \space inputs \space are \space squared $ 1 Normalized NSE where the inputs to the original NSE equation are squared before normalizing.
Log-transformed NSE (NSELog) $ NSELog = 1 - \frac{\sum (\log(Q_o) - \log(Q_s))^2 } {\sum (\log(Q_o) - \log(\bar{Q_o})^2)} $ 1 Log-tranformed Nash-Sutcliffe Efficiency
Weighted NSE (NSEWt) $ NSEWt = \frac{NSE + NSELog}{2} $ 1 Capture flow timing and magnitude errors jointly via the NSE metric and somewhat reduce the peak flow emphasis of NSE by including the log-transformed metric.
Pearson Correlation (Cor) $ r = \frac{ \sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y}) }{%\sqrt{\sum_{i=1}^{n}(x_i-\bar{x})^2}\sqrt{\sum_{i=1}^{n}(y_i-\bar{y})^2}} $ 1 Flow timing.
Root Mean Squared Error (RMSE) $ RMSE = \sqrt{\frac{1}{n}\Sigma_{i=1}^{n}{\Big(\frac{d_i -f_i}{\sigma_i}\Big)^2}} $) 0 Goodness of Fit between observations and model output
Percent Bias (Bias) $ Bias = \frac{\sum{(sim - obs)}} {\sum{obs}} * 100 $ 0 Goodness of Fit between observations and model
Kling-Gupta Efficiency (KGE) $KGE = 1 - \sqrt{ (s[1]*(r-1))^2 +(s[2]*(\alpha-1))^2 + (s[3]*(\beta-1))^2 }$where:$\alpha=\sigma_s/\sigma_o$ $\beta=\mu_s/\mu_o$ 1 Single metric combining timing and magnitude errors
Multi-Scale Objective Function (MSOF) $ MSOF = \sqrt{\sum{(\frac{\sigma_1}{\sigma_k})^2}\sum{(q_{o,k,i} -q_{s,k,l}}(X))^2 } $ 0 Simultaneously consider contributions from a wide range of time scales of aggregation during the calibration process (i.e., mimicking manual calibration), and reduce the likelihood of the search getting stuck in small ‘pits’, by smoothing the objective function surface
hyperResMultiObj $ hRMO = w_o (1.0 - NNSE) + (w_1 \mid P_e \mid ) + (w_2 \mid V_e \mid) \\ \\ where: \\ NNSE = Normalized NSE,\\ P_e = Peas Discharge Error, \\ V_e = Volume Discharge $ 0 Maximize NNSE while minimizing flow volume and peak bias error
peak_bias $ Peak\_bias = \frac{\overline{\mid P_m - P_o \mid}}{P_o}*100) \\ Where: \\ p_m = model peaks for matched events \\ P_o = observation peaks for matched events $ 0 Quantifies the overall peak bias in matched events between the model and observation throughout the calibration window.
event_volume_bias $ Event \space Volume \space bias = \overline{ \frac{\mid V_m - V_o) \mid} {V_o}} * 100 \\ where: \\ V_m = volume \space model \\ V_o = volume \space observed $ 0 Quantifies the overall volume bias in matched events between the model and observation throughout the calibration window.
eventmultiobj $ W1 * abs(Pbiaspeak) + W2 * abs(VEFlow) $ 0 Account for peak error and volume error.
peak_tm_err_hr $ peak\_tm\_err\_hr = \overline{(\mid T_m - T_o \mid)} \\ where: \\ T_m = model \space time \space (hours) \\ T_o = observation \space time \space (hrs)$ 0 Quantifies the central value of time offset between observed and modeled peak events in hours
Probability of Detection (POD) $$POD = \frac{a}{a+c}$$ 1 Probability that a flood was observed when it was forecasted; where 1 is a perfect score and 0 is the worst.
False Alarm Ratio (FAR) $$ FAR = \frac{b}{a+b}$$ 0 Probability that a flood was forecasted, but not observed; where 0 is perfect and 1 is the worst
Critical Success Index (CSI) $$ CSI = \frac{a}{a+b+c}$$ 0 The proportion of correctly forecast floods over all floods, either forecast or observed
Stedinger’s (1981) lognormal estimator of correlation (corr1) $ r1 = \frac{e^{\sigma^{2UV}}-1}{\sqrt{e^{\sigma^{2U}}e^{\sigma^{2V}}}} $ 1 Stedinger’s (1981) lognormal estimator of correlation (corr1) Purpose: Improved estimator of correlation when hydrologic data is skewed (non-normal), which is often the case for daily and subdaily streamflow https://doi.org/10.1080/02626667.2019.1686639
Lamontagne-Barber Efficiency Mixture Model (lbem) See Ref. 1 Improved theoretical estimator of efficiency based on NSE for when hydrologic data is skewed (non-normal) and exhibits periodicity (e.g., seasonality), which is often the case for daily and subdaily streamflow.
Lamontagne-Barber Efficiency Estimators (lbemprime) See Ref. 1 Improved theoretical estimator of efficiency based on KGE for when hydrologic data is skewed (non-normal) and exhibits periodicity (e.g., seasonality), which is often the case for daily and subdaily streamflow. (Lamontagne 2020 https://doi.org/10.1029/2020WR027101)

In the above table, POD, FAR, and CSI are using the threshold predefined in the obsStrData.Rdata file for the threshold.

After 1st iteration model (WRF-Hydro) run is complete, then a script called calib_workflow will do the followings:

Let's take a look at few of these plots.

As seen above some of the plots are duplicated with the tag name outlier. This is added to the calibration workflow to remove the very large numbers (ouliers) from the plots so they are readable. Figures with the tag _outlier are containing all the interations and are the ones that we will check here. Let's begin with checking how the objective function is progressing with iterations.

The above plot keep tracks of how objective function defined in the setup.parm file is progressing with iterations. It should be noted that DDS is a minimization algorithm. If user select an objective function like KGE that the ideal value is the max value, then 1-KGE is used as objective function. The same goes for different variant of NSE, correlation coefficient and lbem. A healthy calibration procedure would asymptote at the higher iterations numbers when we are getting close to the user defined number of iterations. The red star sign shows the iterations which has the best results so far.

Next, we will take a look at how other error metrics are progressing with iterations next.

The above figure shows how different streamflow metrics are evolving with iterations. Ideally the overall trend is that with the improvement of the objective function, all other metrics improve as well. However, that is not the case always and therefore we keep an eye on the performance of the model to make sure we are not degrading other model performance aspects during calibration.

Plot below is the equivalent plot from official NWMv21 calibration. Note that the objective function in NWMv21 was 1 - weighted NSE and log NSE while the objective function used in this exercise is 1 - KGE. As part of NWMv30 RnD, we tested several different objective functions and decided to use KGE for NWMv30 onboarding, and therefore used in this training. Also note, there is a number of metrics that are reported in NWMv30 that did not exists in NWMv21 calibration. In NWMV21 for this gage with improvements in the objective function,other metrics such correlation coefficient, KGE and etc were also improved. This is the desired outcome.

Let's take a look at how parameters have evolved with calibration progress.

The above plot shows how parameters are changing with iterations, it is one way to make sure you are calibrating all the parameters that you flagged. Note that in the early iterations, a larger subset of parameters are being perturbed and as we reach to the end of calibration (larger iterations), fewer parameters are purturbed. It might be more obvious in plot below which is for the same gage from NWMv21.

In initial iterations, the DDS algorithm searches globally and as the procedure approaches the maximum user-defined number of iterations, the search transitions from a global to a local search. This transition from a global to local search is achieved by dynamically and probabilistically reducing the search dimension which is the subset of the calibration parameters that will be updated in a given iteration.

The probability of a parameter to be chosen for inclusion in the search is equal to P(i) = 1 - ln(i)/ln(m), where i is the iteration number and m is the maximum iteration number. Therefore the possibility of a parameter to be chosen reduces with increase in iteration numbers. In the initial iterations, almost all the parameters will be modified and as it approaches the maximum number of iterations it will only modify a few parameters or only one. Parameters selected in each iteration are perturbed within the defined parameter range. The suggested lower and upper limits were shown in lesson 1. The limits are selected based on previous literature review and experts opinion.

The maximum number of iterations used in the previous versions of NWM calibration was set to 300 except for the domains that were too large (> 5000 km2) in that case 150 iterations were used for the calibration.

Next, we will take a look at streamflow hydrograph and the scatter plots.

This plot has only the model simulations with the default parameters (first iteration), the best parameters and the last iterations. We do not plot all iterations in order to have a cleaner picture. However, all the time series are saved in a Rdataset (called proj_data.Rdata) and could be used afterward to plot any other iterations if required. Below plots are the equivalent plots from NWMv21 for the same gage.

Validation Diagnostics:

After calibration is finished, the model runs for both the default and best parameters (from calibration step) for the full duration of (calibration/validation). Note, the validation period does not have to be right after calibration, and in fact validation period could be an earlier time than calibration period. Full Period is the extent that has both calibration and validation period in it. If there is a gap between the two period, the full period will have that gap in it also.

Error metrics are calculated for calibration, validation and full period. Finally a set of diagnostic plots are generated depending on what options user has selected in the setup.parm file. We have only calibrated using streamflow observations in the example provided in lesson 2 and 3, and therefore only streamflow plots are provided here. Let us take a look at the metrics first.

This plot summarizes how model performance has changed from default to best calibration iterations as well as how they perform during an independent time period (validation). Usually the performance of the validation period is not as good as the calibration period.

The main task after calibration and validation workflow finishes is the classification of the calibration basin to Donor, Keep and Drop basins that will be used then in regionalization. Definition of Donor, Keep and Drop basins is as follows:

Below is the same plot from NWMv21. As one can see all the metrics improved compared to the default for both the calibration and validation period.

Let's also look at the hydrograph and scatter plots for our current experiment.

Here is the equivalent plots from NWMv21.

Working with the DataBase

All the paramters and error metrics for each iteraion of calibration and validation are stored in the data base and could be accessed using the sqlite library in python. Since our example case has only one gage, it is very easy to understand the content of the DB. Let's take a look at all the tables that exists in the DB and their content. Let's begin with the Domain_Meta table.

We provided definitions of all the inputs in Lesson 2, refer to the provided table for descriptions. gage_id, link_id and domain_path were the only required inputs from the user, and the rest of the information can be left empty and some are filled out by the python workflow. Our example domain is the contributing area to USGS gage 01447720 and the ComID of link (where gage is located on) is 4185779. The link_id or ComID of reach is identifier used to collect streamflow data from the model outputs. domainID is the identifier that is used to connect this table to other tables in the database. There is a unique domainID for each gage. domaiID is assinged by the python workflow and not the user.

Now let us take a look at the Job_Meta table.

The JobID is specified by the user at the time of submission of the calibration and bundles all the basins that have the same calibration properties. The rest of the information in this table comes from the setup.parm file that is filled by user (refer to Lesson 1). For description of all the columns refer to Lesson 2.

Now let ue take a look at the Job_PARAM table.

The above table lists all the parameters used for calibrations with the range and default value to start from. All the basins in a given job will use the same default and range for their parameters. The table was prepared and explained in Lesson 1. Refer to Lesson 1 for definition of these parameters. Next, take a look at the Calib_Stats table which has the calibration related information.

The table above sumamrizes error metrics for each gage (specified with the domianID) at every iterations during calibration. Parameters for each iterations are stored in a different table called Calib_Params, let us check it out.

The User can define which iterations are best performing from the Calib_Stats and find the correponding parameter values in Calib_Params table.

Finally the last relevant table is Valid_Stats which summarizes the error metrics at the calibration, validation and full period for model runs with the default and best parameters.

Simulations states whether the model run is with the default or calibrated (best) parameters. evalPeriod indicate what the period of evaluations is. The options are calib which is the calibration period, valid is validation period and full is the period (continuous) that has both calibration and validation in it. Note the claibration and validation period do not have to be following each other, there could be a gap between the two or validation could happen at the years prior to calibration.

Visualize Spatial Output Fields

Here we will spatially visualize the difference in a parameter field between the default CTRL output and the BEST selected optimum parameter.

Let us first take a look at how values for LKSATFAC in the fine resolution file (Fuldom.nc) file has changed during our experiment.

As shown above LKSATFAC starts from the uniform value of 1000, this is the default value specified in the calib.parm file. Since this parameter is uniform in space, the parameter is a substitute and not a multiplier, therefore the best LKSATFAC would be a uniform value as well which is equal to what you got from the calibration process.

Next we will take a look at the parameter SMCMAX which exists in both soil_properties.nc file and also HYDRO_TBL_2D.nv file. SMCMAX is a function of soil type and therefore it is not uniform in space (check the CTRL run) and is being formulated as a multiplier in the calibration procedure.

As shown above, the CTRL SMCMAX and BEST SMCMAX are both spatially varying, and the spatial pattern of the CTRL is preserved in the BEST since we used a multiplier. Multiplier value in third figure is the same as the value in the best parameter table above. This parameter is used by both LSM and hydro part of the model and therefore exists in the Hydro 2d files also. Let's confirm they match betweent the two files (they should).

Conclusion:

We reviewed the plots that are generated during calibration and after validation and described the type of the information they provide. Database content also was investigated. You can now move to the next lesson that describes some of the errors that you might run into while trying to use PyWrfHydroCalib tool.